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Amna Noreen and Kare Olaussen 



o 

(N 
c3 



43 

Oh 
i 

43 

-i— > 



(N 
(N 

O 

(N 

> 

X 

c3 



Abstract — The Frobenius method can be used to represent 
solutions of ordinary differential equations by (generalized) 
power series. It is useful to have prior knowledge of the 
coefficients of this series. In this contribution we demonstrate 
that the magnitude of the coefficients can be predicted to 
surprisingly high accuracy by a Legendre transformation of 
WKB approximated solutions to the differential equations. 

Index Terms — Second order ODEs, Regular singular points, 
Frobenius method, Legendre transformation, WKB approxima- 
tion. 



I. Introduction 

RECENTLY we have developed and used code for 
solving ordinary Frobenius type differential equations 
to very high precision (T], 0, 0, like finding the lowest 
eigenvalue of 



ip"(x)+x A ip(x) = sip(x) 



(I) 



to one million decimals. The general class of equations 
treated in Q is of the type 
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Following the Frobenius method our solution is represented 
by a convergent series 



m-\-v 



(3) 



m=0 



where the coefficients a m is generated recursively in parallel 
with a brute force summation of the series. The individual 
terms in ^ may grow very big, leading to huge cancellations 
and large roundoff errors. It is therefore useful to have some 
prior knowledge of the magnitude of the a n 's before a high- 
precision evaluation — to set the computational precision 
required for a desired accuracy of the final result, and to 
estimate the time required to complete the computation. 

We have found that \a m \ can be estimated surprisingly 
accurate from a WKB approximation of the solution, fol- 
lowed by a Legendre transform. For the general class of 
equations Q the WKB integrals and the Legendre transform 
must be done by (ordinary precision) numerical methods. 

In the remainder of this paper we first derive a Legendre 
transform relation between the magnitudes \a m \ and IV'(^)I) 
slightly generalized to take into account a logarithmic cor- 
rection, and next use the WKB approximation to estimate 

^(4 

Manuscript received Mars 18, 2012. 

A. Noreen is with the Department of Physics, NTNU, N-7048 Trondheim, 
Norway, e-mail: Amna.Noreen@ntnu.no. 

K. Olaussen is with the Department of Physics, NTNU, N-7048 Trond- 
heim, Norway, e-mail: Kare.01aussen@ntnu.no. 



II. Legendre transform method of solution 

Our method of solution is based on the hypothesis that 
the sum ^ for large \z\ receives its main contribution from 
a relatively small range of m-values. Introduce quantities u 
and s(m) so that 



_ e s (™)^ 



Our assumption is that 



(4) 



with the main contribution to the sum coming from a small 
range of m-values around a maximum value m. The latter 
is defined so that s'(m) + u — 0, s"(m) < 0. Now write 
m — fh + Am, and approximate the sum Q over Am by a 
gaussian integral. This gives 



In summary, we have found the relations 

u = — s'(to), 



S(u) = s(m) — {y + m)s (m) + — log 
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(6) 



-s"(m). 

This is essentially a Legendre transformation between s(m) 
and S(u). Consider a small change u — > u + Su. To 
maintain the maximum condition we must also make a 
small change m — > m + 8m, with 8m — —8u/s"(m). 
I.e. s"(m) = —u'(m). This is consistent with the result of 
taking the m-derivative of equation ([5j. One further finds 
that So(u) becomes 

1 



S Q (u) + S' (u) 8u + ^(u) 8u 2 + 



Sq(u + Su 

= s(m) + (y + m)u + (m + v) 8u — 



Su 2 



1 

2s"(m 

giving the relations 

(m + v) = S' (u), 

s(m) = Sq(u) — uS' (u), 
s"{m) = S>>(u)-\ 

Equation ^ just says that (dm/du) = (du/dm)^ 1 . We are 
only able to compute S(u) directly, not Sq(u). However, they 
only differ by a logarithmic term, hence we will approximate 

log(-s"(m)) = - log ^'(u) -iogS"'(u). This gives 



(7) 
(8) 
(9) 



So(«)w5(«)-ilog(27rS"(«)) 



(10) 



which can be used in equations (|7j|9]l when we have com- 
puted S(u). 



III. WKB APPROXIMATION 

It remains to find S(u). Here we will use the leading order 
WKB approximation to find a sufficiently accurate estimate. 
When z = is an ordinary point, i.e. when v- — 0, v + = 1, 
the leading order WKB solution to (j2j) is 

y/Qo/Q{x)exp(^J Q(t)dtj, (11) 

where Q 2 (z) — Yl n =i v n zTl ^ 1 > ar >d Qo — Q(0). This repre- 
sents a superposition of the solutions tp±(z). The difference 
between the v + and v_ solutions is at worst comparable to 
accuracy of our approximation; hence we will not distinguish 
between them. 

When z = is a regular singular point we use the 
Langer corrected WKB approximation to obtain leading 
order solutions in the form 



Numerical and estimated coefficients 
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Here Q 2 (z) = \s 2 ( V+ - v_) 2 + £^ = o v ™ z " +i > and Q = 
Q(0). In equation (12i we distinguish between the v + - and 



-solutions, because the difference v + — v_ may in principle 
be large. 

The WKB integrals must in general be done numerically, 
sometimes along curves in the complex plane. This requires 
careful attention to branch cuts. Here we will only give 
some examples where most of the calculations can be done 
analytically. 

A. Example 1: Anharmonic oscillators 
Consider the equation 

d 2 



dy 



:^(y) + (y 2 + c 2 y^(y) = 1 



(13) 



for real c so that c 2 > 0. For large y the typical solution 
behaves like 

- e^ 3+c \ (14) 

neglecting the slowly varying prefactor. For a given value of 
\y\ this is maximum along the positive real axis. Hence, with 
x = y 2 = e u , we find as a leading approximation 

S(u) = | (J u + 3c 2 e^ . 

In this case the Frobenius series can be written 

oo oo 

*(y) = «™ V 2m+V = E A ™(y)> (15) 

m— Tci— 

with v = 0, 1. Ignoring the log(S"('u))-term in (|7] [8] [Tojl 
we find 

m = i (ei" + c 2 e3"V (16) 
log (|a m |) = (| - |«) et" + c 2 (1 - \u) es". (17) 

For c = an explicit representation is 

2 

log |a m | = -m(l - log 2m) . (18) 

This is plotted as the lower curve in figure [Tj It fits 
satisfactory with the high-precision coefficients generated 
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Fig. 1. Comparison of nu merical coef ficients a m (points) with estimates 
(full-drawn lines) based on j!6||17| and l |28||29| . The estimates of log \ a m \ 
are accurate up to corrections which depend logarithmically on m. 



numerically, but there remains a connection which depends 
logarithmically on m. For nonzero c the parametric represen- 
tation provides equally good results, as shown by the upper 
curve in figure [TJ 

The conclusion of this example is that for a fixed (large) 
x we expect the largest term of the power series to be 



max|A m (x)| - e 



(19) 



neglecting a slowly varying prefactor. Further, the maximum 
should occur at 



i (x^ 2 + cV/ 2 ) 



(20) 



Finally, estimates like equation (18 1 for the coefficients a m 
may be used to predict how many terms A4 we must sum to 
evaluate ip(x) to a given precision P, based on the stopping 
criterium 

IflAtl^ < icr p (21) 

As can be seen in figure [2] the agreement with the actual 
number of terms used by our evaluation routine is good, 
in particular for high precision P. But keep in mind that a 
logarithmic scale makes it easier for a comparison to look 
good. 

Next consider the logarithmic corrections. Including the 



prefactor of equation (14i changes S(u) by an amount 

1 



AS(u) 



log (e u 



(22) 



Including the log(5"(u))-term in the relation between S(u) 
and Sq{u) changes Sq by an additional amount 



AS (u) 



ilogQe§" + Jc 2 e^). (23) 



For c 2 = this changes the relation ( 1 8 i to 

log \am\ = I (2m + 5/2) (1 - log (2m + 5/2)) . (24) 



For \a m \ this essentially conesponds to a factor m 
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Predicted and calculated summation lengths 
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Fig. 2. This figure compares the a priori prediction, based on equation 118} , 
of the number of terms .M which must be summed in order to evaluate 
^(y) for c = to a desired precision P with the actual number of terms 
computed by. 



Ratio of numerical to estimated coefficients 
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Fig. 3. This figure shows the ratio between the computed coefficients 
a m and the crude pr edic tion j I s| (labelled |a„'|) and the logarithmically 
corrected prediction 124) (labelled |t4^|)' F° r eas y comparison we have 
in both cases adjusted an overall constant such that the ratio is unity for 
rn = 3. 



B. Example 2: Double well oscillators 

The same procedure also work for the equation 

d 2 



dy 



(25) 



which however is a little more challenging since the maxi- 
mum value of ^(ye'^l sometimes occur for tp ^ 0, i.e. for 
complex arguments. 

For large y the typical solution behaves like 



■ e 3 ' 



(26) 



neglecting the slowly varying prefactor. Equation ( |25[ ) can be 
transformed to the form pb by introducing x — y\^(y) = 
ip{x). Hence, with x = y — e" 

S(u) = max |Re (el^+W _ 3c 2 e^ u+i ^ . 



1 /2 

The maximum occurs for cos \ip = — i(l + c 2 e~ u ) 
when e u > |c 2 , and for cos |<£> = — 1 otherwise. This gives 



S(u) 



c 2 e u/2 

|(e" + c 2 ) 3 / 2 



l e 3u/2 for e u < l c 2 ) 



for e 11 > ±c 2 . 



This implies that 



for e u < |c 2 , 
for e u > |c 2 , 



log(|ohm|) 



l e n/2 ( c 2 _ e „) 
l e « (e^ + c 2 ) 1 / 2 

c 2 e«/ 2 - (| - §«) 



3 "/ 2 for e" < 



(27) 



(28) 

±c 2 
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for e" > |c 2 . 



(29) 

This representation compares fairly well with the numeri- 
cally generated coefficients, as shown by the middle curve in 
figure [T] However, in this case the coefficients a m have a lo- 
cal oscillating behaviour. The representation ( |2"8| [29| ) should 
be interpreted as the local amplitude of this oscillation. 

The conclusion of this example is that we expect the 
largest term of the power series to be term of the series 
to be 

max \A m {x)\ ~ e ^*+^ 3/ \ (30) 

rn 

neglecting the slowly varying prefactor. Further, the maxi- 
mum should occur at 

„2\i/ 2 
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(31) 



IV. Conclusion 

As illustrated in this contribution the coefficients of Frobe- 
nius series can be predicted to surprisingly high accuracy 
by use of Legendre transformations and lowest order WKB 
approximations. We have also tested the validity of the 
method on many other cases. 
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